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Second-order  intensity  statistics  of  a  coherent  signal  in  the  presence  of 
a  random  background 

R.  BARAKATt 

The  statistics  of  a  harmonic  signal  (coherent  component)  mixed  with  a  random 
background  (incoherent  component)  of  a  specified  spectra  profile  (power  spec¬ 
trum)  is  still  a  problem  of  interest.  The  purpose  of  the  present  paper  is  to  study  the 
second-order  intensity  statistics  of  such  a  signal/background  situation  using  the 
generalized  Karhunen-Loeve  expansion  developed  for  use  in  photon  counting  and 
laser  speckle. 


1.  Introduction 

In  a  classic  paper  written  several  years  after  the  actual  work  (done  during  World 
War  II),  Kac  and  Siegert  (1947)  determined  the  exact  first-order  statistics  (i.e. 
probability  density  and  moments)  of  a  square  law  detector  for  a  gaussian  random 
field,  and  for  an  harmonic  signal  buried  in  a  gaussian  random  field.  Their  approach 
involves  the  construction  of  a  homogeneous  integral  equation  whose  kernel  is  the 
covariance  function  of  the  random  field  using  what  is  now  termed  the  Karhunen- 
Loeve  expansion  (Selin  1965,  Thomas  1964);  although  we  can  also  term  the 
construction  the  Kac-Siegert  expansion  in  as  much  as  they  derived  it  independently. 
The  eigenvalues  determine  the  probability  density  function  of  the  detacted  intensity 
(square  of  the  field  amplitude),  while  both  eigenvalues  and  eigenfunctions  are  needed 
for  the  probability  density  of  the  intensity  of  the  signal  and  field.  Emerson  (1953) 
also  discussed  the  problem  from  an  alternative  viewpoint,  using  the  method  of 
cumulants  to  avoid  solving  the  associated  homogeneous  integral  of  Kac  and  Siegert. 
For  further  work  on  the  problem,  see  Slepian  (1958).  As  Mayer  and  Middleton 
(1954)  have  pointed  out,  the  original  expansion  method  of  Kac-Siegert  is  not 
general  enough  to  handle  higher-order  statistics  such  as  product  moments.  However, 
Kac-Siegert  also  presented  another  method  of  solution,  the  ‘direct'  method,  which  is 
capable  of  handling  higher-order  statistics..  The  direct  method  requires  an  appro¬ 
priate  transformation  to  express  the  output  in  terms  of  the  input;  the  statistics  of  the 
output  are  then  determined  by  suitable  additional  transformations  w’ith  respect  to 
the  original  input  statistics,  Mayer’  and  Middleton  exploited  this  approach  to 
determine  the  higher-order  statistics  of  the  output  due  to  a  square  law  detector  for 
both  narrow-band  and  broad-band  inputs.  Kac-Siegert  only  considered  the  broad¬ 
band  situation. 

The  purpose  of  the  present  paper  is  to  restudy  the  above  problems  for  the 
second-order  intensity  statistics  using  a  generalization  of  the  original  Kac-Siegert 
expansion  approach.  I  employ  two  point  detectors  to  interrogate  the  random  input 
(with  and  without  the  harmonic  signal  present).  These  detectors  operate  for  the  same 
time  interval  but  are  delayed  relative  to  each  other  by  a  variable  time.  The  analysis  is 
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performed  via  an  expansion  of  the  random  field  over  these  two  disjoint  time 
intervals,  using  a  generalization  of  the  Karhunen-Loeve  series  developed  for  similar 
problems  in  photon  counting  (Jakeman  1970,  Blake  and  Barakat  1973,  Barakat  and 
Blake  1980)  and  laser  speckle  (Barakat  and  Blake  1978).  In  this  approach,  a 
homogeneous  integral  equation  is  constructed  over  the  two  disjoint  time  intervals 
wherein  the  two  detectors  operate.  The  eigenvalues  and  eigenfunctions  (which  obey 
an  unusual  orthogonality  condition)  are  evaluated  and  used  to  fabricate  a  double 
generating  function;  from  this  double  generating  function  the  various  product 
moments  of  the  integrated  intensities  can  be  obtained  by  differentiation.  Analysis  is 
confined  to  the  narrow-band  situation,  although  there  is  no  difficulty  in  studying  the 
broad-band  situation.  I  feel  that  the  approach  via  the  generalized  Karhunen-Loeve 
expansion  offers  a  more  satisfying  physical  picture  than  the  direct  method  in  as 
much  as  the  detector  time  intervals  and  delays  are  an  inherent  part  of  the  analysis 
via  the  associated  disjoint  integral  equation. 

2.  Formal  solution 

The  complex-value  U{t)  of  the  total  disturbance  is  the  sum  of  a  deterministic  (or 
coherent)  component  U^{t)  and  a  random  background  term  V^jt) 

V{t)=V,{t)-rl^{l)  (1) 

The  coherent  component  is  given  by 

{7c(0  =  ?cexp(-icu,/)  (2) 

where  is  a  constant.  The  random  background  L'^ft)  is  taken  to  be  a  zero-mean, 
complex-valued,  spatially  stationary  gaussian  random  process,  i.e. 

U,U)=U^m  +  iU'^i,)  (3) 


where  L'(;>(/)  is  the  stochastic  Hilbert  transform  of  L"^\t).  Both  and  L'|;>  are 
real-valued.  Now  and  [/{,'*(/)  have  the  same  gaussian  probability  density 

function  (PDF);  also 


<[/[('(/)>  =  <tf(r)>  =  0  (4) 

<f/l,^'(/)O'L'>(/)>  =  0  (6) 

where  erg  is  the  variance  of  U^^{t)  and  —  is  the  corresponding  correlation 
function,  Since  is  a  gaussian  random  process,  it  is  completelv 

characterized  by  its  mean,  variance,  and  correlation  function. 

The  disturbance  L  {t)  is  interrogated  by  two-point  detectors,  the  first  detector 
operating  during  the  time  interval  (“T/2,  7/2)  and  the  second  during  the  time 
interval  (t  — 7/2,  r-f  7/2)  w'here  t  is  the  variable  time  delay.  The  integrated  intensi¬ 
ties  Qj  are 
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“■■i: 


|U.(/)+L/b(0|M/ 


r/2 

T/2 


(7) 


|t/c(0  +  I7bWPd/ 


We  can  expand  Ug[t)  in  a  generalized  Karhunen-Loeve  series  (Jakeman  1970, 
Blake  and  Barakat  1973,  Barakat  and  Blake  1978,  1980)  over  the  disjoint  intervals 


AM-TI2,TI1),  A^  =  {x-TI1,x  +  TI1) 


(8) 


so  that 


t/b(0 


teA^  and  A^ 


0, 


t^A^  and  A^ 


(9) 


The  following  conditions  are  to  be  satisfied. 

(1)  The  {Uy}  are  random  coefficients,  independent  of  /;  and  the  uncorrelated 
gaussian  random  variables  (hence,  independent  random  variables) 

iU,Ur>  =  cld,,  (10) 

where  {ffj}  are,  as  yet,  unknown,  real  non-negative  constants.  It  is  essential 
that  [/b(0  be  gaussian,  for  if  it  is  not  then  the  expansion  coefficients  { will 
not  be  statistically  independent,  although  they  will  still  be  uncorrelated. 

(2)  The  deterministic  functions  {ij/ft)}  are  to  form  a  complete  orthonormal  set 
over  both  Ai  and 

The  precise  statement  of  the  orthogonality  condition  is  quite  unusual.  Consider 
the  weighted  sum  of  the  integrated  intensities:  Aifii  +  where  A, ,  7.2  are  arbitrary 
real,  non-negative  parameters  appearing  in  the  two-fold  generating  function  of  the 
integrated  background  intensities  Qj*” 

Qb(7.i,7.2)  =  <exp(-/ifi'2‘’')>  (11) 


We  have 

;.,^il  +  ;.2n2=|iL/b(0Pd^-Kj)|{7,(/)|M/  +  (^C,(r)O':(r)dr-^(t)U^(/)C/,(/)d/  (12) 


where 


rTi2 

-^7.2 

J  -  T/2 


•t  +  T’2 

t-T/2 


(13) 


We  now  require 

oW0'/'f(0d/  =  ^w  (14) 

The  first  term  on  the  right-hand  side  of  (12)  can  be  evaluated  using  (9)  so  that 
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Ain''’>+A2n<2'>'=  i  \uf 


k  =  0 


(15) 


The  second  term  yields 

(C|[;,(t)|Mr=2i^/r(;.i+A,) 

The  interaction  integrals  (third  and  fourth  terms)  are 


£  l/t(t)(^,(/)exp(iAOdr 

k  =  0 


(C[/*(t)l/,(0d/  =  ^,  £  Ut 


k  =  0 


(l)%t)Qxp{  —  iAt)dt 


(16) 


(17) 


Here 


A  =  (cOc-aJo)  (18) 

is  the  frequency  offset  describing  the  position  of  the  signal  with  respect  to  the 
maximum  of  the  background  power  spectrum.  Consequently  (12)  reduces  to 

k  =  0 


+ I  ,  /2) + c:  z  u,Gt(X, , ;.,)  (19) 

k^O  k=0 


where 


GkO-i ,  Az)  =  C)(;f»i(0exp  (iA/)  dr 


(20) 


To  evaluate  the  unknown  constants  in  (10),  we  must  construct  an  integral 
equation  whose  kernel  is  the  correlation  function  of  U^{t).  The  integral  equation  is 


[a, 

•7/2  r 

+  ^2 

\  J 

-7/2 

r^+Tii 


t-T/2 


(21) 


The  correlation  function  —  of  the  random  background  field  is  given  by 


rbUi  -  h)  =  exp  {(iwo(^  -  h)]9i‘i  -  h)  (--) 

where  ^b(^i“”^2)  is  the  correlation  function  centred  at  zero  frequency  and  cJq  is  the 
frequency  at  which  the  power  spectrum  (lineshape)  of  the  background  radiation  is  a 
maximum.  If  we  set 


(23) 


then  (15)  reads 


/ 

•7,2  r 

+  ^2 

\  J 

-7/2  J 

‘t  +  T/2\ 


-7/2  / 


(24) 
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independent  of  cdq.  This  is  the  basic  integral  equation  for  determining  {ajcf-,  it  is 
not  of  the  standard  Fredholm  type  because  of  the  presence  of  two  disjoint  domains 
of  integration.  A  second  difficulty  is  that  the  eigenvalues  ajja  are  implicit  functions 
of  Aj  and  Aj. 

The  double  generating  function  of  the  random  background,  can  be 

expressed  as  an  infinite  product.  From  (11)  (with  [/.(O^O  for  the  moment)  we  have 

Gb(Ai.A,)=  T-  n  (25) 

Jo  ,  Jo  »=o 

where  d^C/jt  =  d[/j."’dl/i‘’,  and  fl{U^}]  is  the  joint  probability  density  function  of  the 
statistically  independent  gaussian  random  variables 

f[{Uk}]=f\:^^^P{-\Uk\^l<^k}  (26) 

k=0^°k 

Upon  substituting  (15)  into  the  integrand,  we  encounter  a  standard  gaussian 
integral;  the  final  result  is 

Gb(2i,A2)=  n  [l+<r?(Ai,A2)]’'  (27) 

k  =  0 

Now  for  the  double  generating  function  with  the  coherent  signal  included.  Given 
(11)  and  (19),  it  follows  that 

e(Ai,A2)  =  exp{-|c%£'7(A,  +  A2)}  f[  Qki^.u^-i)  (28) 

k^O 

where 
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exp{-(l+c7*  ^)Vl-a^V^}(iV^  = 


exp 


4 


4(1  +  a,- 

^  al-bl 
4(1 +<7,-^) 


f  exp{-(  +  cri^)l^)t  ^}cosfctli;diy*  =  — ^^^exp 
*/  ”■  00 

Thus,  QkQ-iJ.-^  reduces  to 

Consequently,  the  double  generating  function  for  the  background/signal  is 


e(;.i,/2)=exp{-i<^/r(;., +;.,)}  fl  (i+<^.-)''  n  ^p 

k=0  k^O 

=  Qc(‘^-1  >  ■^•2)Qb('^-l  ’  '^•2)Qbc(''-l » •^•2) 


(1+ffJ 


(33) 

(34) 

(35) 


(36) 


where  Q^,  ^re  the  generating  functions  for  the  coherent,  background  components 
respectively  and  Qbc  is  the  generating  function  for  the  interaction  between  back¬ 
ground  and  signal. 

The  product  moments  of  the  integrated  intensities  can  be  obtained  by  differentia¬ 
tion  of  Q{ki,?.2) 


<fl?ns>=(-i)'’*‘’ 


^P  +  9 

a;.?  c4 


^(■^•i' ‘^•2)1*1=22=0 


(37) 


This  completes  the  general  solution. 

3.  Product  moments  of  ft,  small  T  regime 

It  is  obvious  that  the  larger  is  T,  the  greater  is  the  smoothing  of  the  field 
amplitudes.  Consequently,  we  want  to  employ  as  small  a  T  as  possible.  Assume  that 
T  is  small  compared  with  the  distance  over  which  the  background  field  correlation 
function  decays  to  its  exp"'  value.  Under  this  condition,  only  the  k=0  terms  in 
and  Qbc  contribute;  all  /:>0  terms  are  essentially  negligible.  It  has  been  shown  that 
(Jakeman  1970,  Blake  and  Barakat  1971,  Barakat  and  Blake  1978.  1980) 


<7o(7-i  ,  7.2)  =  a^T{/.i+/.2)  +  (T^T\\ -f  |g(/,  -/,)|-)7.i;.,  (38) 

The  Go  function  can  be  evaluated  by  the  mean  value  theorem  for  integrals.  The 
result  is 


Co  =  b{k^T  +  A2T)  exp  (il/  ^  --  t.l)  (39) 


where  i  =  <^o(o)- 

Upon  carrying  out  the  necessary  manipulation 


Q{/.x  ,;.,)= (1  -b  (75) "  ‘  exp  { -  |(?/r(;.i  -i- ;.,)}  e.xp 


KcGoiVs] 

,(I+<T5)  j 


(40) 


Second-order  intensity  statistics  of  a  coherent  signal  225 

The  product  moments  of  the  integrated  intensities  were  obtained  by  differentia¬ 
tion  of  I  employed  the  symbol  manipulation  program,  SMP,  using  a  VAX 

computer  to  effect  the  differentiations.  The  first  few  product  moments  are 

=  =  +  (41) 

= <t1i  -h  iffWl'] + led"  (42) 
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5.  PRODUCT  MOMENTS:  BACKGROUND/SIGNAL 

The  product  moments  of  the  integrated  intensities  were  obtained  by  differentiation 
of  Q(Ai,A2),  Eq.  (4.3),  according  to  Eq.  (3.16);  I  employed  the  symbol  manipulation 
program,  SMP,  using  a  VAX  computer  to  effect  the  differentiations.  The  first  few  product 
moments  are: 

(Oi)  =  {fi2)  =  ^"  +  iec|-  (o.i) 

(Q1Q2)  =  +  g'iT )]  +  2(7‘1^c|'  +  l^cl^  (5.2) 


(Q^Q2)  =  (QiQ2^^2(7''[1  +  25-(t)] 

+  l<fc|* [4  -f  ~g'(T )  +  6  +  26  cos  Ar]  +  3a~  4-  |^c|^  (5-3) 


(QjQ^)  —  4cr*[l  +  45"(r)  +  g^(T )] 

+  <^*’ICc|'{4  +  S^'lrdl  +  6  cos  Ar]  +  86} 

+  cr'^|^c|‘‘[2  +  f/'V)  +  46  +  86  cos  At] 

+4^decr+i6r  (5.4) 

Note  that  the  freciuency  offset  A  only  appears  in  (QjQo)  ^i^cl  (f^ifla). 

In  the  homodyne  case  where  the  frequency  of  the  signal  coincides  with  the  maximum 
of  the  power  spectrum  of  the  background,  then  A  vanishes  and  the  product  moments  then 
depend  upon  ^(t)  only. 
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FILON  TRAPEZOIDAL  SCHEMES  FOR  HANKEL  TRANSFORMS 


OF  ORDERS  ZERO  AND  ONE 


Abstract 


Algorithms  for  evaluating  zero  order  and  first  order  Hankel  transforms  using  Filon 
quadrature  philosophy  are  developed  in  the  context  of  a  trapezoidal  approximation  rather 
than  of  a  Simpson’s  rule  approximation  previously  discussed.  Unlike  the  Filon/Simpson 
algorithm  previously  developed,  the  Filon/trapezoida.1  algorithm  tends  to  saturate,  in  that 
increasing  the  number  of  quadrature  points  does  not  materially  increase  the  accuracy. 
Numerical  examples  are  given  and  discussed. 

1.  Introduction 


Previous  papers  [1,2]  were  devoted  to  the  numerical  evaluation  of  Hankel  transforms 
of  orders  zero  and  one. 


H{r)  =  /  h{p)J„{rp)pdp 


(1.1) 


using  Filon  quadrature  pliilosophy.  In  [l],  the  Filon  approach  is  outlined  in  some  detail  for 
Eq  1.1  with  n  =  0.  There  the  slowly  varying  part  of  the  integrand,  h{p),  is  approximated  by 
a  quadratic  function  over  the  basic  quadrate  panel.  For  n  =  1,  see  [2].  It  was  necessary  to 
consider  h{p)  =  ph{p)  as  the  basic  function  to  be  expressed  as  a  quadratic.  As  with  Filon’s 
original  approach  to  Fourier  integrals,  the  errors  incurred  in  Eq.  1.1  are  proportional  to 
the  derivatives  of  h{p)  and  h{p)  themselves  rather  than  to  the  whole  integrand,  hence  are 
relatively  independent  of  r. 

In  many  areas,  we  do  not  require  great  accuracy  for  the  Hankel  transforms,  but  need 
to  maintain  a  given  accuracy  more  or  less  uniformly,  independent  of  the  magnitude  of  r. 
The  purpose  of  the  present  communication  is  to  develop  the  trapezoidal  version  of  the 
Filon-Hankel  approach.  In  [1,2],  the  scheme  is  reall}'’  a  generalization  of  Simpson’s  rule, 
since  we  are  employing  a  quadratic  approximation  of  h[p)  and  h{p)  over  the  panels.  The 
present  approach  is  essentially  a  generaliza.tion  of  the  trapezoidal  quadrature  scheme,  since 
we  are  approximating  h{p)  and  h{p)  by  a  linear  function  over  the  panels. 
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2.  Trapezoidal  Algorithm/Zero-Order  Transform 


Consider  the  integral 

rvt+t 

Hk{r)  =  /  h{p)Jo{rp)pdp 

''Jfk 

for  two  points  pk+i  and  pi,.,  where 


Pt+i  -Pk  =  S 


(2.1) 


The  points  where  k  =  0,1,. 

■  ■  ,N,  are  a  subset  of  [a, 6]. 

Approximate  h{p)  by  a  straight 

line  between  pk  pat+i: 

h{p)  =  A  +  Bp 

(2.2) 

It  follows  that 

A  =  -{pk+ihk  -  Pkhk+i) 

(2.3) 

B  =  ^(/lA+i  -  hk) 

(2.4) 

where  h{pk)  =  hi,. 

4 

Upon  setting  u  =  h{p)  and  dv  =  J^,[Tp)pdp,  we  integrate  i7/:(r)  by  parts;  the  end 
result  is 


Hk{r)  =  -[hk+iJi{rpk+i)pk+i  -  hkJi{rpi,)pk]  -  ^[h'k+M'^Pk+l)  -  /ijtSij(rpi-)l  (2-5) 

where  we  have  used 

J  yJi){y)dy  =  xJi{x)  (2.6) 

The  numerical  evaluation  of  the  function 

S„(i)  =  /  yJi{y)dy  (2.7) 

J\) 

is  discussed  at  some  length  in  the  Appendix  of  [l],  to  which  we  refer. 

Since 

=  K  =  ]{fn+i  -  hk)  (2.8) 

then  Eq.  2.5  reduces  to 

Hk{r)  =  -  [hk+iJi{rpk^i)  -  hkJi{rpk)]  -  —{hk+i  -  hjt)  [Sii(rpji.+i )  -  Sii(rpt)]  (2.9) 
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To  evaluate  the  integral  over  [a,  6],  i.e.,  Eq.  1.1,  divide  [a,  6]  thusly 


b  =  a  + NS 


(2.10) 


Hence 


The  final  result  is 


N 

ffw  = 

k-{) 


(2.11) 


H{r)  =  -[hi^+iJi{Tpr,r)  -  hoJi(rp(,)]  -  -  hk)[$u{rpk+i)  -  SoC^pife)]  (2.12) 


k=[) 


This  is  the  basic  formula  for  the  Filon/trapezoidal  scheme  for  zero-order  Hankel  transforms. 
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3,  Numerical  Example  for  the  Zero-order  Case 


As  with  the  Filon/Simpson  algorithm  for  the  zero-order  Hankel  transform,  we  consider 
as  in  [l] 

=  ^[arcos(p) -p(l 0<p<l  (3.1) 

H{t)  =  [ — 0  <  r  <  oo  (3.2) 

T 

as  our  test  case. 

Unlike  the  Filon/Simpson  algorithm,  the  Filon/trapezoidal  algorithm  tends  to  satu- 
rate,  in  that  increasing  N  does  not  materially  increase  the  accuracy  of  the  quadrature. 
Of  course  this  is  not  surprising,  hecause.  we  are  now  approximating  h[p)  by  straight  lines 
rather  than  by  quadratic.  In  the  context  of  our  numerical  example,  perhaps  the  best  way 
to  see  this  is  to  fix  r  while  varying  examining  the  absolute  error 

l^f(7’).-.\-:..  t  -  F/(7-)o.mimtHcl!  (3.3) 

Some  typical  values  of  the  absolute  error  are  displa3’’ed  in  Table  1.  This  table  does  not 
require  any  detailed  comment. 

In  evaluating  the  various  .T-Bessel  functions,  we  employed  Mason’s  algorithm  [3],  which 
is  probably  the  most  accurate  currently  available. 
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4.  Trapezoidal  Algorithm/First-Order  Transform 


Analagous  to  Eq.  2.1,  we  write 

l‘Pi+1  . 

Hk{r)  =  /  h{p)Ji{Tp)dp 
Jjn- 

where 

h(p)  =  h[p)p 

Now  set  u  =  h(p)  and  dv  =  Ji[rp)dp.  Upon  integrating  by  parts,  we  obtain 


Hk{r)  =  -  -[h{pk-ki)Jn{rpk+i)  -  Hpk)Ji){rpk)] 

iLfj;'/ 

^2 


-r  -^[h'{pk+i)Ri){rpk+i)  -  h' {ph)R„{rpk)] 


where  we  ha.ve  used  the  indefinite  integral 


and 


J  Ji{y)(iy  = 

0 

=  /  Ji){y)dy 
J[) 


See  Appendix  A  of  [2]  for  the  numerical  evaluation  of  this  function. 
As  in  Section  2,  we  can  sum  the  various  panels,  leading  to 


=  -  -  [fhxMrp.x)  -  fill  J(i(rp,i)] 
1  ^ 


(4,1) 


(4.2) 


(4.3) 


(4.4) 


(4.6) 


(4.6) 


This  is  the  basic  formula  for  the  Filon-trapezoidal  scheme  for  first  order  Hankel  transforms. 
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5.  Numerical  Example  for  the  First-order  Case 

Let  us  consider  the  example  [4] 

/i(p)  =  p(l  -  ,  0<p<l  (5.1) 

H{r)  =  0  <  r  <  oo  (5.2) 

2t 

The  numerical  results  are  summarized  in  Table  2.  As  with  the  corresponding  algorithm  in 
Section  2,  this  algorithm  also  tends  to  saturate  as  N  increases. 

6.  Summary 

The  Filon/trapezoidal  scheme  is  not  meant  to  be  a  direct  competitor  to  the 
Filon/Simpson  scheme.  The  main  purpose  of  this  quadrature  scheme  is  to  maintain  a 
given  accuracy  (provided  it  is  not  too  extreme)  more  or  less  uniformly,  independent  of 
the  magnitude  of  the  indejjendeiit  variable.  An  added  advantage  of  this  scheme  is  its 
speed  of  execution,  an  important  aspect  for  such  problems  as  beam  propagation  in  an 
inhomogeneous  or  random  medium,  where  the  integral  must  be  computed  a  large  number 
of  times.  Reference  is  made  to  [5]  for  invention  of  the  Hankel  transforms  using  the  sampling 
expansion  in  connection  w'ith  the  Filon/Simpson  and  the  Filon/trapezoidal  schemes. 
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Table  1.  Absolute  error  of  Filon/trapezoidal  calculations  for  zero-order  case,  example  in  Eqs  3.1  and  3.2. 


r 

N  =  100 

N  =  200 

iV  =  300 

N  =  400 

6 

3.727E-06 

9.202E-07 

4.069E-07 

2.282E-07 

12 

3.063E-06 

7.722E-07 

3.447E-07 

1.945E-07 

18 

2.099E-06 

5.418E-07 

2.443E-07 

1.386E-07 

24 

1.083E-06 

2.928E-07 

1.343E-07 

7.690E-08 

30 

1.901E-07 

6.871E-08 

3.427E-08 

2.048E-08 

40 

8.976E-07 

2.354E-07 

1.068E-07 

6.082E-08 

50 

5.229E-07 

1.199E-07 

5.167E-08 

2.861E-08 

60 

1  2.053E-07 

6.985E-08 

3.411E-08 

2.015E-08 

80 

1  2.544E-07 

4.877E-08 

1.941E-08 

1.026E-08 

Table  2.  Absolute  error  of  Filon/trapezoidal  calculations  for  first-order  case,  example  in  Eqs  5.1  and  5.2. 


r 

N  =  100 

N  =  200 

N  =  300 

A^  =  400 

6 

2.764E-04 

9.789E-05 

5.330E-05 

3.463E-05 

12 

2.186E-04 

7.801E-05 

4.260E-05 

2.772E-05 

18  1 

1.769E-04 

6.422E-05 

3.527E-05 

2.302E-05 

24 

i 

1.360E-04 

5.080E-05 

2.822E-05 

1.851E-05 

30 

9.406E-05 

3.715E-05 

2.096E-05 

1.386E-05 

40 

1.129E-04 

4.214E-05 

2.331E-05 

1.526E-05 

50 

1.029E-04 

3.578E-05 

1.927E-05 

1.243E-05 

60 

7.152E-05 

2.157E-06 

1.093E-05 

6.818E-06 

70 

2.863E-05 

4.011E-05 

9.669E-07 

2.174E-07 

80 

1.479E-05 

1.215E-05 

7.884E-06 

5.537E-06 

10 


NUMERICAL  EVALUATION  OF  FOURIER  INTEGRALS: 
FILON  QUANDRATIVE  VERSUS  THE  FFT 


1.  INTRODUCTION 


We  are  concerned  with  the  numerical  evaluation  of  the  finlite  range  Fourier  integral 


F{p)e'^P 


dp, 


—  '00  <  V  <  00 


(1.1) 


assuming  that  F(p)  varies  slowly  compared  to  the  complex-valued  exponential  term.  The 
case  where  F(p)  is  of  rapid  variation  is  analyzed  by  a  different  approach,  see  Section  5.  For 
large  values  of  |v|,  a  graph  of  the  integrand  consists  of  positive  and  negative  areas  of  nearly 
equal  size.  The  addition  of  these  areas  results  in  substantial  loss  of  accuracy.  Filon  [1] 
conceived  the  idea  of  retaining  Simpson's  rule,  but  requiring  that  only  F(p)  be  fitted  to  a 
quadratic  over  the  basic  subinterval  instead  of  the  entire  integrand  jF’(p)exp(ivp).  The  fact 
that  only  F{p)  has  to  be  approximated  means  that  the  number  of  subintervals  to  be  taken 
can  be  relatively  small  in  many  cases  of  practical  interest.  An  additional  feature  of  the 
Filon  quadrative  philosphy  is  that  the  error  incurred  is  relatively  independent  of  v  because 
the  error  is  proportional  to  the  derivatives  of  F(p)  itself  rather  than  to  F(p)exp(ivp). 


2.  FILON/SIMPSON  ALGORITHM 

We  will  now  sketch  in  some  detail  the  Filon/Simpson  algorithm  for  evaluating  the 
finite  range  Fourier  integral  (under  the  fundamental  assumption  that  F{p)  varies  slowly 
compared  to  the  complex-valued  exponential  terms  also  in  the  integrand. 

Consider  the  integral 


fkiv) 


f 


■F{p)e‘"’'  dp 


(2.1) 


which  is  effectively  a  double  panel  of  three  quadrature  points:  p2k+2-:  'P2k+ii  P2k  where 


^  =  iP2k+2  —  P2k+1 )  =  (P2jt-t-l  —  P2fc)" 


(2.2) 


This  panel  is  a  subset  of  the  larger  interval  [a,  6]. 

Following  Filon  we  assume  that  F{p)  is  smooth  enough  to  be  approximated  by  a 
quadratic  function  over  the  interval  pok  <  p  <  P2fc-i-2 

F{p)  =  -f  boip  -  P2k+i )  +  ?>3(p  -  P2k+i )'  (2.3) 

where  the  b's  are  as  yet  unknown.  Upon  solving  for  the  b's  we  have 


^1  =  Fok 

h  =  ■^{F2k+2  -  Fojfc) 

^3  =  r^{F2k+2  ~  ~F2k+l  +  F2k)  (2-4) 

where  F2k+2  =  F{p2k+2)i  etc.  The  first  and  second  derivatives  of  F{p)  are  also  needed. 
Differentiating  Eq.  (2.3)  and  then  employing  Eq.  (2.4)  jdelds 

^2k+2  —  l^{^F2k+2  —  4F2fc+i  +  F2k)  (2-5) 

^2k  —  '^i~F2k+2  +  4F2fc+i  —  3F2k) 

Finally 

^2k  =  -p{F2k+2  -  2F2fc+i  +  F2k)  (2.6) 
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Next  integrate  Eq.  (1.1)  twice  by  parts,  the  final  result  is 


fk(v)  =  -  ^2^+26*"^=*+^) 

V 

V“  -  I  - 

+  4F4(e‘''”‘+“  -  e'*'”*) 


(2.7) 


Substitution  of  Eqs.  (2.5)  and  (2.6)  ultimately  leads  to 


fk{v)  =  — (E2fc+2e*''P^^+^  -  ^2^6*"^^-) 

V 

+  ne  ,(3E2Jfc+2  ~  4E2Jt+l  +  F2k)F''^^-’'~-  +  {F2k+2  ~  4E2jt+l  +  3E2Jfe)e”’^^’‘ 

2ov- 

+  -^^F2k+2  -  2E2jfe+a  +  F2ife)(e‘"^'=''+=  -  )  (2.8) 

This  expression  is  the  basic  building  block  of  the  Filon-Simpson  algorithm. 

To  evaluate  Eq.  (1.1).  we  divide  the  interval  of  integration  [a,  6]  thusly 


h  =  a  +  NS 


(2.9) 


where  N  is  an  even  integer  (as  in  the  standard  Simpson  method;  thus 

(N-2)/2 

f{v)  = 

fc=0 

Combining  terms  in  the  summation,  this  eventually  becomes 

(N-2]/2  , 


fiv)  =  -{FoF^P°  -  Fnc^^^-  +  V  ^  +  i 

V  ^  \  28v-  6-v^ 


fc=0 


where 


.10) 


(2.11) 


Go  =  F2  —  4Ei  +  .3Eo 

Gfc  =  F2k+2  —  4F2fc+i  +  6E2fc  —  4E2fc_i  +  F2k-2 

Gjv-2  =  3Ejv  —  4Ejv-i  +  Ff^-o 
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(2.12a) 

{2.12b) 

(2.12c) 


Ho  =  —F2  +  2jFi  —  Fq 


(2.13a) 

(2.13b) 

(2.13c) 


ajid 


Hk  =  —F2k+2  +  ~F2k  +  l  —  -F2k-1  +  F2k 


Hn-2  =  -  2i^iV-l  +  Fn-2 


In  Eqs.  (2.13b)  and  (2.13c),  we  have  0  <  k  <  (N  —  2)/2.  Equation  (2.11)  is  the  basic 
expression  for  the  Filon/Simpson  quadrature  algorithm. 

A  modification  must  be  made  when  u  =  0,  because  of  the  presence  of  v  in  the  denom¬ 
inators  of  Eq.  (2.11).  Fortunately  this  is  straightforward  because 


m 


- 

O  a 


Fl,p)  dp. 


(2.14) 


Since  F{p)  is  assumed  to  be  fairly  smooth,  the  integral  can  be  evaluated  directly  using  a 
standard  quadrature  formula;  our  case.  Simpson's  rule. 

Equation  (2.11)  constitutes  the  corrected  complex  exponential  version  of  the  Filon 
scheme  found  in  [2].  Equation  (2.11)  is  markedly  different  than  the  final  forms  of  the 
sine  and  cosine  versions  derived  by  Filon  [1].  He  combined  terms  in  a  manner  that  was 
conducive  to  hand  calculation  techniques  of  his  time,  resulting  in  calculations  that  were 
independent  of  F{p)  and  so  could  be  tabulated  without  the  knowledge  of  F{p).  However, 
computational  abilities  are  considerably  advanced  now,  so  that  it  is  unnecessary  to  store 
tables  of  values  in  order  to  perform  the  integration.  Direct  calculation  using  Eq.  (2.11)  by 
computer  is  not  difficult.  Reference  is  made  to  [3.4]  for  the  explicit  e.xpressions  of  the  sine 
and  cosine  versions  of  Filon’s  original  analysis:  they  different  significantly  in  form  from 
Eq.  (2.11)  particularly  with  respect  to  the  very  troublesome  a  and  j3  terms  in  the  original 
version. 

We  omit  any  discussion  of  the  error  incurred  as  it  is  still  a  matter  of  contention  [5]. 
However,  the  error  is  proportional  to  the  derivatives  of  F{p)  itself  rather  than  to  the  entire 
integrand,  hence  are  relatively  independent  of  the  magnitude  of  the  variable  v. 
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3.  FILON/TRAPEZOIDAL  ALGORITHM 

In  some  situations  we  do  not  require  great  accuracy  but  need  to  maintain  a  given 
moderate  accuracy,  more  or  less  uniformly,  independent  of  the  magnitude  of  v.  To  this 
end  we  develop  a  trapezoidal  version  of  the  previous  algorithm,  the  Filon/ trapezoidal 
algorithm. 

We  now  sketch  this  algorithm.  Consider  two  points  pk+i  and  pk  which  are  a  subset 
of  [a,  6].  Approximate  F(p)  as  a  straight  line  between  pk+i  and  pk 

F{p)  =  A  +  Bp,  Pk<p<Pk+i  (3.1) 


Here 


where  S  =  (pk+i  -  Pk)-  Also 


^  =  -^{Pk+lFk  —  PkFkJrl 

B  =  “  ^k) 


F'  =  “  ^k)- 


Upon  integrating  by  parts  and  using  the  alDove  equations,  we  have 


fk(v)  =  -  -  Fke^'-’^’') 


(3.2) 


(3.3) 


(3.4) 


To  evaluate  the  integral  over  [a,  6]  i.e..  Eq.  (1.1).  divide  [a,  6]  thusly  b  =  a  +  N6  so 


that 


The  final  result  is 


N 


/(»■)  = 


(3.5) 


k=0 


/(v)  =  -  -  Fue'^^’o) 

V  ' 

1  ^ 


(3.6) 


fc=0 


When  V  =  0,  we  again  employ  Eq.  (2.14) 
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We  should  emphasize  that  the  Filon/trapezoidal  algorithm  is  not  meant  to  be  a  direct 
competitor  to  the  Filon/Simpson  algorithm  as  regards  great  accuracy.  Rather  its  main  use 
is  to  maintain  a  given,  but  moderate  accuracy,  independent  of  the  magnitude  of  when 
speed  of  execution  is  important. 

In  the  special  case  where  the  limits  of  integration  of  Eq.  (1.1)  are  symmetric  (i.e., 
a  =  —6),  then  the  Filon/trapezoidal  algorithm  can  be  written 

/(»)  =  (f)  E 


where  (2M  +  1)  is  now  the  number  of  quadrature  points,  and 


